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ABSTRACT 


The  Institute  for  Naval  Oceanography  (INO)  has  established  the  Experimental  Center 
for  Mesoscale  Ocean  Prediction  (ECMOP)  to  facilitate  research  and  development  of  ocean 
models.  ECMOP  aims  to  take  advantage  of  existing  infrastructures  to  provide  support  to  Navy 
and  academic  researchers  in  an  end-to-end  model  evaluation  effort.  In  addition,  if  and  when 
necessary,  ECMOP  performs  internal  research  and  development  as  an  adjunct  to  this  support. 
ECMOP  is  composed  of  three  functional  subsystems,  or  modules:  the  Verification  Module 
(VERMOD),  the  Visualization  Module  (VISMOD),  and  the  Data  Sub-System  (DASS).  This 
technical  memorandum  describes  the  current  state  of  VERMOD,  in  particular,  its  statistical  and 
physical  criteria  for  model  verification,  and  the  software  structural  characteristics  of  the  module. 
Finally,  we  discuss  several  of  the  verification  concepts  anticipated  as  future  enhancements  of 
VERMOD. 


ACKNOWLEDGMENTS 

We  would  like  to  thank  Mr.  James  Corbin  for  his  comments  on  an  earlier  draft  and 
Ms.  Lydia  Harper  for  her  editorial  help. 


1.  INTRODUCTION 


The  Institute  for  Naval  Oceanography  (INO)  has  established  the  Experimental  Center  for 
Mesoscale  Ocean  Prediction  (ECMOP)  to  provide  capabilities  that: 

•  permit  development,  demonstration,  and  evaluation  of  mesoscale  ocean  prediction  models; 

•  furnish  information  pertinent  to  the  transition  of  ocean  modeling  systems  to  Navy  sponsors 
(and  "decision  makers")  in  an  objective,  orderly,  and  efficient  manner;  and 

•  are  modular  in  nature  and  can  be  transitioned  as  appropriate  to  operational  Navy. 

While  adhering  to  operational  concepts  developed  by  the  Naval  Oceanography  Command, 
these  capabilities  are  being  realized  through  a  modular  structure  that  emphasizes  a  standard¬ 
ized  interface  among  data,  numerical  models,  and  processing  systems.  ECMOP  is  comprised 
of  three  functional  sub-systems:  the  Verification  Module  (VERMOD)  provides  objective  eval¬ 
uation  routines,  the  Visualization  Module  (VISMOD)  defines  user  definable  graphic  routines, 
and  the  Data  Sub-System  (DASS)  provides  data  acquisition  and  management.  In  all  areas  of 
ocean  modeling,  whether  the  effort  is  R<SiD  or  operational,  there  is  a  need  for  basic  technical 
support  capabilities  for  common  functions  to  acquire,  manage  and  analyze  data;  in  situ,  re¬ 
motely  sensed,  or  model  generated.  All  these  functions  incorporate  some  aspects  of  three  basic 
technical  capabilities:  (1)  data  base  management;  (2)  graphical  and/or  visualization  schemes; 
and  objective  evaluation  routines.  The  role  of  ECMOP  is  to  focus  on  these  three  basic  tech¬ 
nical  capabilities  commensurate  with  the  needs  of  the  R&D  community  by  way  of  building 
them  into  a  state-of-the-art  facility  using  leading  edge  technology  in  software  and  hardware. 
At  the  same  time,  ECMOP  is  to  evaluate  these  state-of-the-art  capabilities  for  applicability 
to  the  operational  community  and  when  identified,  make  them  available  for  transition.  The 
conceptual  details  of  ECMOP  were  initially  given  by  Leese  (1988, a, b).  The  present  ECMOP 
configuration,  in  its  modular  form,  is  shown  in  Fig.  1,  wherein  the  three  modules  communicate 
externally  using  the  network  common  data  form  (netCDF)  jackets. 

This  technical  memorandum  documents  the  current  status  of  VERMOD.  VERMOD  was 
originally  developed  with  a  minimal  number  of  statistical  and  physical  tools  for  model  verifi¬ 
cation  and  evaluation  in  order  to  substantiate  the  software  design  characteristics.  These  tools 
are  accessed  through  an  intuitive  graphical  user  interface  (GUI)  that  requires  little  user  indoc¬ 
trination  or  training.  The  GUI  functions  within  the  X-windows  environment  as  implemented 
under  the  UNIX  operating  system.  Its  features  are  accessible  to  heterogeneous  environments 
and  across  computer  networks  and  platforms.  Enhancements  and  modifications  will  be  made 
to  VERMOD  according  to  user  requirements.  Guidance  for  crucial  enhancements  has  been 
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derived  from  the  INO  Summer  Colloquium  on  Model  Evaluation  (INO,  1989),  as  well  as  input 
from  the  community  of  scientific  users. 

The  statistical  measures  included  in  VERMOD  1.0  are  the  Root  Mean  Square  (RMS) 
error  with  both  systematic  and  unsystematic  decompositions  (Wilmot  et  al.,  1985),  pattern 
correlation,  and  two  criteria  for  gauging  the  departure  from  normality,  skewness  and  kurtosis. 
These  techniques  are  described  in  the  text  and  illustrated  with  computer  simulated  examples. 

Section  2  gives  a  brief  description  of  the  software  architecture.  Description  of  the  statistical 
measures  is  presented  in  Section  3,  which  is  followed,  in  Section  4,  by  three  examples  that 
illustrate  the  concepts.  Finally.  Section  5  briefly  examines  possible  future  enhancements  to 
VERMOD. 

2.  VERMOD  ARCHITECTURE 

The  VERMOD  application  embodies  characteristics  drawn  from  the  other  modules  within 
ECMOP.  The  software,  its  design  and  the  supporting  hardware  are  briefly  discussed  below. 

2.1  Software 

All  ECMOP  software,  including  VERMOD,  is  designed  to  execute  under  the  UNIX  oper¬ 
ating  system.  To  function  properly,  VERMOD  must  utilize  services  provided  by  the  remaining 
ECMOP  modules.  Therefore,  an  understanding  of  VERMOD  requires  a  description  of  ECMOP, 
both  in  its  present  configuration  and  that  configuration  to  which  it  is  evolving  as  a  complete 
software  package. 

The  GUI,  which  controls  all  modules,  is  implemented  within  X-windows.  The  core  of 
DASS,  the  data  support  module,  is  ’’Empress",  a  proprietary  relational  database  application 
that  implements  the  structured  query  language  (SQL)  and  features  a  capability  to  handle  bulk 
data.  Access  to  the  data  is  via  the  Navy  Environmental  Operational  Nowcast  System  (NEONS), 
an  interface  to  "Empress”  written  in  C  to  support  high  level  access  to  the  data  stored  within  the 
data  base.  NEONS  is  described  by  Computer  Sciences  Corporation  (1991).  The  VISMOD  is 
currently  built  upon  the  National  Center  for  Atmospheric  Research  (NCAR)  application  NCAR 
Graphics.  Model  verification  and  evaluation  within  the  VERMOD  are  performed  using  estab¬ 
lished  and  widely  accepted  quantitative  and  qualitative  evaluation  techniques  as  implemented 
through  FORTRAN  subroutines  and  C  functions.  As  a  complete  package,  ECMOP  software  is 
written  in  "C"  and  "FORTRAN”  with  the  support  of  SQL  and  C-shell  scripts.  The  X-windows 
GUI  uses  the  Open  Software  Foundation’s  ’’Motif’  ("OSF/Motif’)  widget  set.  These  fea¬ 
tures  give  the  ECMOP  software  increased  transparency  across  networks  and  portability  across 
a  variety  of  platforms  that  support  X-windows. 
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2.2  Hardware 


The  ECMOP  software  has  been  developed  on  Sun  Microsystems  hardware.  It  is  easily 
portable  to  other  Unix-compatible  platforms  such  as  those  manufactured  by  Hewlett-Packard, 
Silicon  Graphics  and  Cray.  This  portability  is  made  possible  because  "Empress”,  IMCAR  graphics 
and  X-windows  are  available  on  these  and  other  platforms. 

2.3  Conceptualization 

The  ECMOP  design  is  evolving  into  a  modular  structure.  Currently,  ECMOP  modules  are 
loosely  linked;  however,  redesign  efforts  will  require  user  control  through  the  GUI,  which  oper¬ 
ates  as  the  front  end  to  the  ECMOP  System  Control  Module  (SCM).  The  rationale  underlying 
ECMOP  redesign  is  the  need  to  support  interaction  among  the  several  modules.  The  ECMOP 
SCM  will  integrate  functional  control  of  the  modules  while  enabling  communications  among 
them.  All  existing  modules  are  to  be  enhanced.  In  addition,  a  User  Help  Module  (UHM)  will 
provide  context  sensitive  assistance  to  the  user  at  any  point  within  the  application.  The  DASS 
is  being  redesigned  into  a  Data  Management  Module  (DAMM)  that  will  include  such  features 
as  on-line  directories,  conversion  utilities  and  NEONS  data  ingestion.  A  Model  Control  Module 
(MCM)  is  to  be  added  which  will  allow  ocean  models  to  execute  within  the  ECMOP  envi¬ 
ronment.  The  VISMOD  will  be  upgraded  to  include  capabilities  that  allow  three-dimensional 
visual  analysis.  Enhancements  to  VERMOD  are  discussed  in  this  document. 

3.  STATISTICAL  MEASURES 

In  this  section  we  present  the  derivations  of  the  statistical  measures  that  are  included  in 
VERMOD  1.0,  viz,  the  decomposition  of  root  mean  square  error  (RMSE)  into  systematic  and 
unsystematic  components,  the  pattern  correlation  and  the  two  simple  measures  of  departure 
from  normal  distribution. 

3.1  Preliminaries 

We  would  like  to  compare  two  fields,  m  and  y.  To  compare  the  model  output  with 
observed  data,  we  may  denote  the  model  output  by  the  symbol  m  and  the  observational  data 
by  the  symbol  y.  For  a  proper  comparison,  the  two  fields  must  correspond  to  each  other  in 
some  respect.  This  correspondence  may  be  in  terms  of  a  time  series  for  fixed  spatial  coordinates 
so  that  one  may  compare  over  time  to  obtain  a  spatial  comparison.  Another  simple  example 
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would  be  to  compare  model  output  fields  at  two  different  points  in  time.  In  this  case  the  two 
symbols  represent  the  same  field,  i.e.,  the  model  output,  and  the  comparison  yields  a  temporal 
difference.  The  more  common  need  will  be  to  compare  model  output  with  observed  data. 
This  requires  that  the  two  fields  correspond  in  spatial  distribution  to  provide  a  model-data 
comparison  at  a  given  time.  However,  this  required  correspondence  between  model  output  and 
data  usually  does  not  exist,  since  the  data  are  not  usually  available  at  the  model  grid.  Some 
numerical  manipulation  may  be  necessary  to  achieve  this  correspondence.  When  the  resolution 
of  observational  data  is  finer  than  the  model  grid,  as  in  the  case  of  Multi-Channel  Sea  Surface 
Temperature  (MCSST),  some  (optimum)  interpolation  of  the  data  may  be  performed  to  align 
the  two  fields  to  the  model  grid.  Otherwise,  which  is  usually  the  case,  the  model  output 
data  are  interpolated  to  observation  locations  to  enable  a  valid  comparison.  This  realignment 
introduces  an  error  in  addition  to  those  present  in  either  the  model  output  or  cne  observational 
data. 

Note:  VERMOD  1.0  provides  capability  of  RMSE  computations  when  the  two  fielded  variables 
(model  output  and  observations)  are  already  properly  aligned  to  the  same  grid  frame. 

We  use  the  following  notation.  The  mean  values  of  m  and  y  are  denoted  by  fim  and 
Hy  such  that  S{m)  =  Hm  and  S{y)  =  fiy,  where  £  is  the  statistical  expectation  operator. 
Corresponding  variances  are  defined  by 

crl,  =:  £{m  -  nmf  and  =  £{y  -  ^y)'^ .  (3.1) 


3.2  The  RMS  Criterion 

The  RMS  error  criterion  is  easily  understood  because  its  square,  the  mean  square  error 
(M5E),  is  associated  with  the  variance  and  the  RMSE  provides  a  direct  measure  of  the  difference 
between  the  two  fields  being  compared.  We  first  define  the  MSE  between  two  variables,  m 
and  y,  to  be 

MSE  =  £{m  —  ijY  =  [  {m  —  yf  f{uj)du)  (3.2) 

Ja 

where  /(cu)  is  the  probability  density  function,  or  a  weighting  function  defined  on  the  space 
{u;  €  over  which  m{u>)  and  y{uj)  are  defined.  For  discrete  grid  matched  variables  m  and 
y,  we  can  approximate  the  mean  square  error  as 

.V 

MSE  ~  —  y,)~  (3.3o) 

1=1 
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where  Wi  are  the  weights  such  that  Weights  can  be  varied  to  emphasize 

subsets  within  regions  or  interest.  For  example,  in  the  North  Atlantic  basin  one  may  choose  to 
emphasize  the  Gulf  Stream  by  specifying  larger  weights  for  that  region.  For  equal  weighting 

Wt  = 

N 

MSE  (3.36) 

t=i 


3.2.1  Systematic  and  unsystematic  RMSE 

Decomposition  into  systematic  and  unsystematic  RMS  errors,  as  discussed  by  Wilmot  et  al. 
(1985),  has  been  quite  popular  in  oceanography,  and  was  recommended  for  implementation  in 
VERMCD  by  the  Working  Group  on  Statistical  Procedures  during  the  INO  Summer  Colloquium 
on  Model  Evaluation  (1989).  The  systematic  error  component  identifies  any  bias  between  the 
two  fields  being  compared,  and  the  unsystematic  component  measures  the  goodness  of  fit. 

The  motivation  for  this  decomposition  derives  from  two  common  comparison  scenarios: 
(a)  the  model  output  and  the  data  are  comparable  within  some  random  error,  and  (b)  the 
model  output  is  within  a  fixed  bias  plus  random  error.  For  the  case  (a)  we  will  write 

rrii  -  yi  +  Cj 

where  is  the  random  error.  For  the  case  (b)  we  can  express  the  relationship 

rni  =  00  +  yi  +  e. 

where  0o  's  a  fixed  bias.  These  two  conditions  can  be  generalized  to  the  following  linear 
relationship; 

THi  =  00  +  0iy  +  e,.  (3.4) 

From  this  the  RMSE  decomposition  is  obtained  as  follows.  Let  m  be  the  ordinary  least  squares 
estimate  of  m  computed  from  the  regression  of  m  on  y.  Then,  the  coefficients  0j,  j  =  1,2 
are  estimated  by  minimizing  the  sum  of  squares  —  0o  —  0\yY ■  Let  ho  and  h\  be  the 

minimizing  estimates  of  these  coefficients.  The  predicted  value  is  given  by: 

mi  =  6o+6iy,.  (3.5) 

The  sum  of  squares  of  the  differences  between  the  two  fields  can  then  be  written  as: 

N  y 

-  yif  =  ^iu,((mi  -  m,)  +  (m,  -  y,)]^ 

j=i  1=1 

y  y  y 

-  ^u;,(m,  -  mif  +  ^u’i(m,  -  y,)’*^  -I-  -  m,)(m,  -  y,). 

i=i  i=i  1=1 
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As  shown  in  the  Appendix,  the  cross-product  term  ~  —  Hi)  =  0-  Thus, 

the  above  equation  can  be  written  as  a  sum  of  two  components,  the  systematic  and  the 
unsystematic  components,  given  by: 

;V  N 

-  ijif  =  -  m,f  (rn,  -  y,)T  (3.6) 

1=1  .=  1 


The  unsystematic  component  is  given  by 


RMSE„ 


■  s 

.1=1 


I 


and  the  systematic  component  is  given  by 


RMSE, 


/V 


.«=! 


(3,7) 


(3.8) 


Note  that 

RMSE^  =  RMSE^  -h  RMSE^.  (3.9) 

Wilmot  et  al.  define  RMSE^  as  the  linear  bias  and  RMSE„  as  a  precision  criterion.  The 
above  decomposition  can  be  performed  at  any  depth,  Z .  Suppose  the  model  output  data  m  is 
available  in  the  a  coordinate  system  while  the  other  field,  y,  is  available  at  the  Z  level.  Here, 
the  model  data  at  each  gridpoint  are  interpolated  in  the  vertical  at  the  Z  level  using  a  spline 
function,  thus  leading  to  grid-to-grid  RMSE  evaluation  using  (3.6). 


Note:  In  their  formulation,  Wilmot  et  al.  assume  there  is  no  measurement  error  in  y.  Without 
this  assumption,  the  least  squares  estimation  leads  to  a  complicated  statistical  analysis,  and 
the  final  solution  is  not  that  clean. 


.y.2.2  Grid-to-observations  (option  currently  not  available) 

In  this  case  the  observations  yi,t  =  l,...,iV  at  a  particular  Z  level  are  not  coincident 
with  the  model  grid.  To  compute  RMSE  the  model  output  are  first  interpolated  in  the  vertical 
at  the  Z  level.  A  further  interpolation  is  then  performed  to  derive  model  output  values  at  the 
observation  points  prior  to  RMSE  calculations. 
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3.3  Pattern  Correlation 


Pattern  correlation  is  used  to  estimate  how  closely  the  pattern  of  one  field  resembles 
that  of  another.  Again,  if  and  yi,i  =  are  the  values  of  the  two  fields,  then  the 

coefficient  of  pattern  correlation  is  given  by 

=  -  m)(y,  -  y) 

This  concept  tacitly  identifies  if  there  is  a  linear  relationship  between  the  two  fields.  A  high 
positive  value  of  p  indicates  that  the  patterns  of  highs  and  lows  of  the  two  fields  match  and 
that  the  two  fields  are  almost  linear  transformations  of  each  other.  A  high  negative  value,  on 
the  other  hand,  indicates  that  the  highs  of  the  one  field  correspond  to  the  lows  of  the  other. 


The  pattern  correlation  measure  is  connected  with  the  Wilmot  type  decomposition  by  the 
following  relationship: 


N 

RMSE„  =  —  ihif' 

i-\ 


<T?(i  -  p"). 


(3.11) 


Equation  (3.11)  says  that  the  model  prediction  m  can  be  described  by  m  more  precisely  if  the 
pattern  correlation  is  large. 


3.4  Departure  from  Normality  Measures 


During  the  Trial  Ocean  Prediction  Experiment  (TROPE),  Waters  et  al.  (1990)  recom¬ 
mended  to  include  in  VERMOD  measures  of  skewness  and  kurtosis  that  provide  information 
on  how  closely  the  differences  between  two  fields  resemble  normality.  Skewness  is  defined  as 


and  kurtosis  is  defined  as 


where  pj  is  the  yth  central  moment  computed  from  the  field  differences.  Skewness  tells  whether 
the  differences  are  asymmetrical;  a  positive/negative  value  indicates  that  the  distribution  is 
skewed  to  the  right/left.  The  Kurtosis  measure  indicates  the  flatness  of  the  distribution.  For 
a  normal  distribution  71  =  0  indicating  symmetrical  distribution,  and  79  =  1.  A  negative 
value  of  72  shows  that  the  distribution  is  more  peaked  than  the  normal  distribution,  i.e.,  the 
distribution  is  concentrated  more  around  the  smaller  differences. 
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4.  EXAMPLES 


To  illustrate  the  use  of  the  RMSE  decomposition  and  pattern  correlation  we  will  use  sea 
surface  height  data  from  a  nowcast  for  the  Gulf  of  Mexico  as  the  model  output.  The  observation 
data  will  be  simulated  using  three  sets  of  and  /3i  by 


Vi  =  (mi  -  0o)/f^i 

followed  by  generation  of  the  model  output  by  using  (3.4).  The  noise  added  in  (3.4)  is  Gaussian 
with  zero  mean  and  standard  deviation  <7„=  1  cm.  The  three  sets  of  (/So,/?i)  will  provide 
interpretation  for  three  different  situations,  one  corresponding  to  zero  systematic  bias  and  two 
situations  corresponding  to  when  it  is  non-zero. 

4.1  The  Three  Cases 
Case  A;  RMSEj=0 

This  situation  arises  when,  in  (3.4),  0o  —  0  and  =  1.  For  the  ideal  estimation  situation, 
6o  0  and  6i  «  1.  Then  the  predicted  model  values  are  given  by 


ihi  w  yi. 

Thus; 

N  N 

RMSE^  =  «  ^w,(y,  -  y,f  =  0. 

1=1  j=i 

The  unsystematic  component  in  this  case  is  given  by; 

N  yv 

RMSE^  =  -  naif  ^  ^Wi(m,  -  y^f . 

1=1  i=l 

Simulations  for  this  case  were  performed  using  /?o  =  0  and  /9i  =  1.  The  estimated  values 
of  these  parameters  (Table  1)  are  6o  =  —0.008  and  hi  —  0.996,  which  are  close  to  the  actual 
parameter  values.  The  simulated  model  output  and  observations  are  presented  in  Figs.  2a  and 
2b.  The  sea  surface  height  patterns  in  the  two  are  almost  identical  but  for  the  fact  that,  due 
to  the  simulation  effect,  the  model  output  appears  ragged. 
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Table  1:  Results  of  Simulation  Examples 


{00- 

(/'<), />i) 

(<7„,<T,,) 

Total  RMSE 

P 

Skewness 

Kurtosis 

A 

(0.00, 

1.00) 

(-0.001,  1.000) 

(0.010,  0.001) 

0.011 

0.997 

-0.016 

3.457 

Bl 

(0.10, 

1.00) 

(0.099,  1.000) 

(0.010,  0.099) 

0.100 

0.997 

-0.011 

3.430 

B2 

(0.20, 

0.80) 

(0.199,  0.800) 

(0.010,  0.252) 

0.252 

0.997 

-2.272 

9.352 

The  computed  values  of  RMSEu  is  0  010,  the  same  as  cr„,  while  RMSE^  is  0.001,  which 
is  quite  close  to  zero.  The  RMSE  decomposition  is  in  keeping  with  our  expectation. 

The  computed  value  for  pattern  correlation  in  this  case  is  fairly  high,  0.997.  The  magnitude 
of  the  correlation  is  a  function  of  the  noise  level  (standard  deviation).  Since  the  noise  level  is 
the  same  for  the  three  cases,  the  computed  pattern  correlation  has  a  constant  value. 

Note  that  the  RMSE^  is  0  only  if  both  conditions,  0o  —  ^  arid  I3\  =  1,  are  satisfied.  To 
see  this  we  note  that  if  /?o  =  0  but  ^  1,  then 

mi  «  ^iVi, 

which  leads  to  RMSE^  =  (/?i  -  1)^23, ^  0- 

On  the  other  hand,  if  /3o  0  and  /3i  =  1,  then  m,  ss  /3o  +  J/i  which  yields 

N  N  N 

RMSE^  =  +  y,  -  Vif  =  ^0. 

i=l  i=l  j=l 


Case  B:  RMSE^  ^  0 

The  non-zero  systematic  RMSE  component  is  marked  with  either  /3o  ^  0  or  ^  1. 
When  /?o  7^  0.  the  case  of  interest  is  =  1.  All  other  cases  of  RMSE^  0  are  covered  under 
3\  ^  1.  These  two  situations  are  described  below. 

Case  Bl:  RMSE.  4^  0  with  /ji  =  1 

This  requires  0.  This  is  the  pure  bias  situation  wherein  the  model  output  has  a  fixed 
bias,  /fo,  and  a  random  error,  Pj.  In  this  case  (3.4)  is  written  as: 

m,  -  00  +  y,  +  f,. 


1  2 


(3.14) 


This  case  was  discussed  above  and  yields  RMSE^  ~  =  '^^o-  Thus,  RMSE^  ~  3o, 

as  one  would  expect.  The  simulations  for  this  casa  were  performed  using  So  —  0.10  and 
S\  —  1.  The  computed  values  for  their  estimates  are  ho  =  0.099  and  6]  =  1.00.  As  indicated 
earlier,  the  computed  value  of  the  pattern  correlation  is  0.997,  showing  a  very  good  fit  between 
the  simulated  model  output  and  the  observations;  again,  as  expected.  The  simulated  model 
output  and  data  are  shown  in  Figuies  3a  and  3b.  The  two  show  similar  patterns.  However, 
a  close  examination  reveals  that  the  contour  values  are  displaced  by  Sq  =  0.10.  Note  that 
RMSEu  =  <7,i,  while  RMSE.,  =  0.099  ~  So  =  1- 

Case  B2:  RMSE„  0  with  Ji  4^  1 

In  this  case  Jo  nnay  or  may  not  be  zero.  We  simulated  this  case  with  Jo  =  0.20  and 
Ji  =  0.80.  The  computed  estimates  come  out  to  be  bo  =  0.199  and  bi  =  0.80.  Again, 
we  expect  the  RMSE,,  to  be  close  to  cr„.  In  this  case  the  computed  value  matches  exactly. 
However,  since 

.\  .V 

RMSE^  =  -  ij,y  =  +  (^i  - 

1=1  1=1 

depends  on  both  bo  and  bi,  unlike  the  previous  two  cases,  it  is  difficult  to  accurately  predict 
the  computed  value  of  RMSEj  when  Ji  ^  1.  In  our  simulations,  its  computed  value  is  found 
to  be  0.252. 

Sea  surface  height  patterns  for  the  simulated  model  output  and  observations  are  presented 
in  Figs.  4a  and  4b.  Even  though  the  two  patterns  are  linearly  related,  it  is  difficult  to  discern 
the  relationship  from  such  a  graphical  presentation.  The  RMSE  decomposition  into  systematic 
and  unsystematic  components,  along  with  analysis  of  the  regression  coefficients,  bo  and  bi,  is 
required  to  obtain  the  necessary  insight. 

4.2  Skewness  and  Kurtosis 

Skewness  and  kurtosis  measures  are  computed  for  the  differences  d,  =  m,  —  y,.  Whereas 
the  kurtosis  measu  e  is  difficult  to  interpret,  the  skewness  can  be  quite  revealing.  Skewness 
provides  a  measure  of  symmetry  of  the  differences  in  the  distribution,  d,.  Note  that  this 
distribution  is  symmetrical  only  if  Jj  =  1,  as  in  the  case 


m,  —  Jo  +  J/i  +  ff 


1  3 
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Fig.  4a:  Sea  Surface  Height  contours  ^or  simulated  model  output;  Case  B2:  .fo  =  0.20.  =  0.80 
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This  covers  the  two  scenarios  of  interest:  (1)  when  there  is  no  bias,  i.e.,  /3o  =  0,  and  (2)  when 
there  is  a  fixed  bias,  do-  This  is  evident  from  Figures  5a  and  5b  which  present  the  histograms 
for  the  simulation  cases  A  and  Bl.  When  /?1  ^  1,  we  have 

di  =  do  -f  idi  -  '^)yi  +  e,. 

Here,  the  distribution  of  d,  is  centered  at  however,  its  symmetry  is  determined  by  the  sign 
of  (/?i  —  1).  The  distribution  is  skewed  to  the  left  if  the  sign  is  negative  and  to  the  right 
when  positive.  In  Case  B2,  the  sign  is  negative  leading  to  a  distribution  skewed  to  the  left  as 
confirmed  by  the  histogram  in  Figure  5c.  This  is  also  apparent  from  the  computed  values  of 
skewness  in  Table  1. 

5.  SOME  FUTURE  ENHANCEMENTS 

The  next  version  of  VERMOD  will  incorporate  several  enhancements,  including  a  test 
of  hypothesis  setup  for  the  RMSE  decomposition,  space-time  spectral  decomposition,  and 
empirical  orthogonal  functions  (EOF’s). 

We  will  provide  the  capability  to  assign  confidence  bars  to  the  linear  coefficients  do  and 
di  of  RMSE  decomposition.  This  will  lead  to  an  objective  procedure  of  judging  the  systematic 
and  unsystematic  components. 

We  will  also  develop/assemble  a  methodology  for  space-time  spectral  decomposition  of 
gridded  fields.  This  will  assist  in  determining  a  model's  agreement  with  observed  ocean  fea¬ 
tures,  component  wise.  Although  spectral  decomposition  does  not  provide  a  direct  model  to 
data  comparison  like  the  RMSE  measure,  the  component-wise  verification  can  increase  the 
validity  of  the  model.  Such  applications  require  identification  of  well-established  phenomeno¬ 
logical  features  of  ocean  dynamics  that  a  model  incorporating  proper  physics  and  forcings  must 
produce.  Examples  of  such  phenomenon  often  occur  in  the  literature,  e  g.,  (1)  26-day  oscilla¬ 
tions  observed  in  the  Western  Indian  Ocean  which  were  reproduced  by  Kindle  and  Thompson 
(1989)  using  monthly  averaged  wind  fields  to  drive  the  ocean,  and  (2)  the  eddy  kinetic  energy 
in  certain  frequency  bands  in  the  North  Atlantic  Gulf  Stream  region  (Schmitz  and  Holland, 
1982). 

Empirical  orthogonal  function  (EOF)  representation  of  the  time  series  of  various  oceano¬ 
graphic  parameter  fields  is  another  useful  procedure  of  reducing  the  dimensionality  and  ex¬ 
tracting  information.  EOF  representation  has  been  used  extensively  in  meteorological  and 
oceanographic  analyses.  Recent  examples  from  the  literature  include  Hurlburt  et  al.  (1990) 
and  Carnes  et  al.  (1991)  that  give  useful  applications  to  oceanography. 
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Fig.  5a;  Histogram  of  (model  output  -  observations)  differences;  Case  A;  /fo  =  0,  0\ 
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5b;  Histogram  of  (model  output  -  observatiorrs)  differences;  Case  Bl:  =  0.10,  /f)  =  1.00 
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Fig.  5c:  Histogram  of  (model  output  -  observations)  differences;  Case  B2:  /^d  ^  0.20,  li\  —  0.80 


5.1  Confidence  Intervals  for  RMSE  Decomposition 


Having  decomposed  the  RMSE  decomposition  into  systematic  and  unsystematic  compo¬ 
nents,  it  is  now  necessary  to  anchor  the  methodology  in  an  objective  sense  by  specifying  the 
p— values  along  with  the  estimated  coefficients  bo  and  b^.  Equivalently,  we  must  provide  confi¬ 
dence  intervals  (error  bars).  In  this  way,  it  is  possible  to  determine  the  degree  of  conformation 
to  one  of  the  four  classes;  (1)  no  bias  Case  A,  (2)  simple  linear  bias  as  in  Case  Bl,  (3)  a  simple 
linear  relationship  in  which  /9i  0,  and  finally,  (4)  wherein  the  situation  is  not  described  by  a 

linear  relationship. 

5.2  Spatial  and  Temporal  Spectral  Analysis 

Spectral  decomposition  of  a  space-time  series  is  an  important  analytical  tool  often  used 
in  analyzing  statistical/dynamical  behavior  of  the  ocean.  A  few  examples  of  the  phenomena 
that  can  be  studied  using  space-time  spectral  analysis  are  (William  Johns,  1989  -  personal 
communication):  meander  scales  of  the  Gulf  Stream,  their  propagation  speeds  and  their  growth 
rates.  However,  we  must  be  cautious  in  applying  it  because  a  space-time  series  is  defined  in  a 
two-dimensional  space.  It  poses  greater  difficulties  than  if  it  were  defined  on  a  one-dimensional 
space.  Space-time  series  data  can  be  analyzed  in  three  different  ways:  (1)  analyze  the  data  at 
a  fixed  point  in  space  as  a  simple  time  series;  (2)  analyze  the  data  at  a  single  time  point  as  a 
series  in  space,  being  cautious  in  this  case  to  define  the  spatial  axis;  and  (3)  consider  the  data 
as  a  joint  space-time  series  whose  analysis  provides  wavenumber-frequency  spectra  from  which 
we  can  analyze  traveling  and  standing  waves. 

A  spectral  decomposition  of  a  space-time  series  a(x,f),  with  x  as  space  index  and  t  as 
the  time  index,  is  described  following  Pratt  (1976)  and  Halliwell  and  Mooers  (1979,  1983). 
Let  a{xi^tj)  be  discrete  samples  at  points  Xi  at  times  tj  where  Xi  and  tj  are  equally  spaced, 
such  that  Xi  =  iAx  and  tj  =  jAt,i  =  1,  j  =  1,...,A^.  These  series  represent  the 

amplitudes  of  the  process  being  observed.  We  can  then  define  three  autocorrelation  functions. 
Two  of  these  corresponding  to  time  and  space,  respectively,  are; 


x„ 

1  ^ 

Rioi^hij)  =  ^(0  t  ){M  ^  ^  (^-2) 


N 


0){N 


-ma(xi ,  tj )a{xi ,  tj  -|-  I  (5.1) 

'j=0 
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These  autocorrelation  functions  are  normalized  respectively  by  the  variances  C(x,,0)  and 
C(0,tj).  The  third  is  the  space-time  autocorrelation  function,  defined  by 


M 


N 


^m')  — 


—  ml  1 


■y^  -ma{x,,tj)a{xi  -|-  +  r^).  (5.3) 


Frequency  autospectra  can  be  computed  by  Fourier  transformation  of  the  temporal  autoco¬ 
variance  function  (5.1).  Assuming  the  spectrum  is  space-invariant,  we  can  average  it  over 
its  spatial  coordinates.  Similarly,  the  wavenumber  autospectra  can  be  computed  by  Fourier 
transforming  the  spatial  autocovariance  function  (5.2).  Averaging  over  the  time  domain  can 
be  performed  using  the  assumption  that  the  wavenumber  spectra  are  time  invariant. 

To  compute  the  wavenumber-frequency  spectra,  the  space-time  series  a(x,.fj)  is  first 
expanded  into  temporal  Fourier  harmonics  as 


N 

„(xj)cos(<T„fj)  -F  S„(xi)sin('T„tj)].  (5.4) 

n=l 

The  wavenumber  autospectra  and  cross  spectra  of  the  functions  C„(x;)  and  5n(xj),  for  the 
frequency  an  and  the  wavenumber  km,  are  given  by  WmiCn),  Wm{Sn),  KmiCn,  Sn)  and 
Qm(C'iv,  S/v),  where  Wm  is  the  the  wavenumber  autospectrum  for  Cm  and  Sm,  and  Km 
and  Qm  are  the  wavenumber  cospectrum  and  quadrature  for  C„  and  Sn-  The  ordinary  wave 
number-frequency  spectrum  is: 

E{±km,CTn)  =  \[Wm{Cn)  +  Wm{Sn)]  ±  ^Qm(C„,  5„),  (5.5) 


and  the  total  wavenumber-frequency  autospectrum,  without  regard  to  the  sign  of  k,  is 

T{km,an)  =  ^[Wm{Cn)  +  Tym(5„)].  (5.6) 

The  contribution  from  propagating  waves  is  defined  to  be  the  difference  between  E{+km,o’n) 
and  E(-km,(rn)  given  by 

P{km,(Tn)  =  QmiCn,Sn). 
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5.3  Empirical  Orthogonal  Functions 

Empirical  orthogonal  functions  (EOFs)  are  often  used  to  reduce  the  dimensionality  of  the 
variables  under  consideration.  For  instance,  consider  the  ocean  model  output  field  {a(f,x)  ; 
t  =  1, ...  ,n;  X  =  1, ...  ,p}  such  that  t  is  the  time  index,  x  is  the  space  (grid  point)  index,  and  p 
is  large.  In  the  Chervin-Schneider  (1976)  framework,  we  want  to  test  the  hypothesis  for  change 
in  the  mean  sea  surface  temperature  for  the  grid  under  consideration.  Thus,  the  hypothesis 
pertains  to  the  entire  grid  and  not  individual  gridpoints.  It  is  well  known  that  testing  point  by 
point  does  not  yield  the  proper  significance  at  which  we  are  testing;  a  few  individual  significance 
cases  is  not  equivalent  to  the  significance  testing  of  the  entire  grid  under  consideration,  and 
vice  versa.  In  such  situations  we  can  take  advantage  of  the  unique  decomposition  of  space-time 
data  as  follows.  Let  a(t)  =  [a(l,f), . . .  ,a(/), t)]',  where  prime  indicates  matrix  transposition, 
and 

n 

A  =  ^a(i)a'(().  (5.7) 

<=1 

Let  fc  =  1, . . .  ,p  be  the  p  eigenvectors  corresponding  to  the  eigenvalues  A^.  Then  a(t,  x) 
can  be  expanded  into  a  modal  time  series 

p 

bk{t)  =  a(t)  =  ^d>jt(x)a(x,  f),  (5.8) 

1=1 

where  </>jt(x)  are  the  elements  of  the  eigenvector  Pk,  and  A*  are  the  variances  of  bk{t).  Usually, 
a  few  Afc  account  for  a  large  fraction  of  the  variance  and  the  remaining  A  it  are  negligible.  Modal 
components,  bk{t),  for  only  nonnegligible  A*  are  retained  and  the  remaining  are  ignored,  thus 
reducing  the  dimensionality.  Preisendorfer  et  al.  (1981)  provide  objective  rules  for  selecting 
the  number  of  EOFs.  Relation  (5.8)  can  be  inverted  to  obtain  the  space-time  decomposition 
of  the  field  a(f,  x)  as: 

p 

a{t,x)  =  'Y^bk{t)(^k{x).  (5.9) 

With  this  reduced  dimensionality  we  can  now  perform  ail  statistical  tests  of  hypotheses  that 
were  possible  with  the  original  space-time  series  a{t,x). 
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6.  CONCLUDING  REMARKS 


The  current  capabilities  of  the  Verification  module  VERMOD  being  developed  in 
ECMOP/INO  have  been  described.  Included  in  VERMOD  are  several  of  the  capabilities  of 
VISMOD,  the  visualization  module.  The  two  modules  combined  together  provide  an  effective 
tool  to  the  scientific  community  for  some  aspects  of  model  evaluation.  Version  1.0  includes 
only  rudimentary  statistical  measures  to  quantify  model-data  differences.  The  interpretation 
and  use  of  these  tools  are  illustrated  by  examples  in  which  model  output  and  observational 
data  are  simulated  using  computer  simulations.  This  current  version  is  being  transitioned  to 
the  Naval  Oceanographic  Office. 

The  modules,  VERMOD  and  VISMOD,  by  their  nature,  will  be  in  a  perpetual  state 
of  enhancement.  As  indicated  in  Section  5,  several  enhancement,  e.g,  space-time  spectral 
decomposition  and  empirical  orthogonal  functions  software,  will  be  incorporated.  The  future 
versions  will  be  developed  as  newer  and  more  sophisticated  softwares  are  included  according 
to  the  scientists’  requirements  and  availability/development  of  the  advanced  techniques.  It 
is  anticipated  that  the  majority  of  future  evaluation  techniques  will  be  developed  within  the 
modeling  groups  themselves,  and  integrated  into  VERMOD  where  applicable. 

APPENDIX 

We  show  that  the  cross-product  term 

'^Wi{mi  -  mi){mi  -  yi)  =  0.  (.41) 

;=1 

Note  that  this  cross-product  term  is  a  sum  of  two  terms: 

N  N 

Ti  +T2  =  'Y^w,{mi  -  rhi)mi  -  '^Wi{mi  -  mi)yi. 
i=l  1=1 

We  will  show  that  each  of  these  two  terms  is  zero.  We  start  by  obtaining  the  estimates  bo  and 
61  of  3o  and  3i  by  minimizing  the  sum  of  the  squares  —  0o  —  Then 


m,  =  60  +  hyi. 


(.42) 


The  minimization  yields  the  following  normal  equations; 

N 

-  bo  -  biiji)  =  0, 

1=1 
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(A3) 


N 


Y^Wi(mi  -  bo  -  biyi)yi 

1=1 

Using  (A2)  we  can  rewrite  the  normal  equations  as; 


=  0. 


N 

^^Wi{Tni  —  TTli)  =  0, 
i=l 

N 

'Y^Wi{mi  -  mt)yi  =  0. 

i=l 


Note  that  (A4)  is  the  same  as 


Ti  =  0. 

Multiply  (A3)  by  bo  and  (A4)  by  bi  and  add  to  obtain: 

N  N 

0  =  ^u>i(mi  -  7hi)bo  +  ^u;i(m,  -  m,)biy, 

i=l  i=l 

N 

-  -  m,)(bo  +  bjy.) 

i=i 
N 

=  ^iui(mj  —  mi)mi 
1=1 


=  32. 


(A4) 

(A5) 


(A5) 
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